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I. INTRODUCTION 

The accurate calculation of time-dependent quantum 
observables in correlated electron systems is crucial to 
achieve further progress in this field of research, where the 
vast majority of the computational efforts in the past have 
mainly focused on time-independent quantities. This issue 
lies at the core of the study of a broad array of physi- 
cal phenomena in transition metal oxides and nanostruc- 
tures such as electronic transport, optical excitations, and 
nonequilibrium dynamics in general. Accurate studies of 
time-dependent properties will advance the fields of spin- 
tronics, low dimensional correlated systems, and possibly 
quantum computing as well. For a list of recent efforts 
on these topics by our group, and concomitant set of refer- 
ences for the benefit of the readers, see [1-5] and references 
therein. 

The purpose of this paper is to present an explicit imple- 
mentation of the time evolution within the density matrix 
renormalization group (DMRG) method [6, 7]. Knowledge 
of the widely discussed DMRG algorithm to compute static 
observables is here assumed. Readers not familiar with 
the method are referred to published reviews [8-11], and 
to the original publications [6, 7] for further information. 
Our main focus here is to provide a detailed description of 
the implementation of the time-step targetting [12] algo- 
rithm, and the discussion of a few applications. We also 
provide full open-source codes and additional documenta- 
tion to use those codes. [43] 



The present work builds upon considerable previous ef- 
forts by other groups. In particular, we will mainly follow 
Manmana et al. in Ref. [13]. The time-step targetting pro- 
cedure was also reviewed in Ref. [14]. Since it would not 
be practical to describe in a short paragraph the consid- 
erable progress achieved in this field of research in recent 
years, the interested reader is encouraged to consult the 
aforementioned reviews along with, e.g., Ref. [15], for a 
historical account of the development of the methods used 
in the present publication. 

Since our aim is to discuss a generic method applicable 
to any Hamiltonian and lattice geometry, here we do not 
discuss or implement the Suzuki- Trotter method [16, 17], 
but focus instead only on the Krylov method [18] for the 
time evolution, as described in Ref. [13]. Because its im- 
plementation can be isolated, the Krylov method can be 
applied in a generic way to most models and geometries 
without changes, which is not the case for other methods, 
such as the Suzuki- Trotter method. The readers interested 
in the Suzuki- Trotter method should consult the ALPS 
project [19], where the time evolving block decimation is 
implemented [16]. 

Our goal is to compute observables of the form 

((/.i|e^^%,,(o)Ai,,(i).--A,_i,,(,_i)e-*^*|02). (1) 

where \(j)i) and \4>2) denote generic quantum many-body 
states. This category of observables is sufficiently broad 
to encompass most time-dependent correlations, as repre- 
sented by a number a of local operators ^o.ir(o) ^i,Tr(i) ' ' ' 
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^a-i,ir(a-i) acting on sites 7r(0), 7r(l), 7r(a — 1) of a finite 
lattice, where 7r(j) denotes the lattice site on which the 
operator ^i^7r(i) acts, and the extra index i indicates that 
the operator can be different at each site. 

An immediate application of this formalism and code is 
the study of the evolution of a system that is brought out 
of equilibrium by a sudden excitation. This sudden exci- 
tation can be simulated by the state of the system given 
by the vectors and |02)- Depending on the problem, 
sometimes it is more convenient to assume that the states 
remain unchanged but that it is the Hamiltonian H(t) that 
changes with time. Another application entails the com- 
putation of time-dependent properties of systems in equi- 
librium, such as the Green's function Gij(t). 

The organization of this paper is the following. Sec- 
tion II explains in detail the Krylov method for time evo- 
lution within the DMRG algorithm, focusing on its imple- 
mentation. Section III A applies the method to the case of 
one-site excitations, showing a simple picture of the accu- 
racy of the method. Section III B extends to two-leg ladder 
geometries the results obtained using tight-binding chains, 
employing holon-doublon excitations for the specific study. 
Performance issues are studied in Sec. IV, while Sec. V 
summarizes our results. The first two appendices contain 
derivations of exact results used in the paper. The last 
appendix explains briefly the use of the code, and points 
to its documentation. 



II. METHOD AND IMPLEMENTATION 
A. Lanczos Computation of the Unitary Evolution 

To carry out the previously described program [44] of 
computing observables of the type given by Eq. (1), the 
first goal is to calculate = exTp{iHt)\(j)). The Lanc- 

zos technique [20] provides a method to tridiagonalize H 
into V'fTV, where T is tridiagonal and V is the matrix of 
Lanczos vectors. If the number of those Lanczos vectors 
is n/, and the Hilbert space for |0) has size n, then T is 
a square matrix of size n/ x n/, and V is, in general, a 
rectangular matrix of size ni x n. 

V'^TV cannot be used everywhere as a substitution for 
H , without inducing large errors. But if we start the 
Lanczos procedure [21] with the vector (instead of us- 
ing a random vector as is most frequently done), then we 
can use that substitution accurately in the multiplication 
However, this is not enough here, because we need 
to compute the exponential of H . For this purpose, it has 
been shown [22-24] that « exp(iri)y [</)) with an 

accuracy that increases as time t decreases for fixed n/. 
We will assume that we have taken t small enough such 



that we can regard the expression above to have the accu- 
racy of the Lanczos technique, which is usually high. In 
other words, if t is small enough, we will assume that us- 
ing eyiY>{iTt)V\(t>) as a replacement for exp(zi?i)|0) is 
not worse than using V^TV\4)) as a replacement for H\(j)). 
This will be enough for our purposes, but for details on 
the scaling and bounds of the errors made in each case as 
a function of t and n;, see, e.g., Refs. [22, 23]. 

Since we started the Lanczos recursive procedure with 
the vector then J^j ^■'jl</')j Sj'^a- Finally, we need 
to diagonalize T = S^DS into a ni x ni diagonal matrix 
D with diagonal elements di'. This last step is not com- 
putationally expensive, since T is a x matrix, as was 
noted before. 

Putting it all together, we arrive to 

m)h= E ylkTL'Slj,e''^''S,jTi,oVoM,^ (2) 

k,k' ,1,1' ,j 

for small times t, where the equal sign should be under- 
stood to be valid within the accuracy of the Lanczos tech- 
nique [24]. How to deal with larger times t will be ex- 
plained in section II C. 



B. Targetting States with the DMRG Algorithm 

It appears that now we could use Eq. (2) to compute 
l^i(t)) from some vector \4>i), and likewise \(j)2{t)) start- 
ing from some vector \4)2)- Then we would just apply the 
operators Ao^7r(o) ^i,7r(i) • • • ^a-i,Tr(a-i) to those states 
within a DMRG procedure, to achieve our aim of com- 
puting Eq. (1). But the DMRG algorithm is not imme- 
diately applicable to arbitrary states, such as \(l)i(t)), and 
was originally developed to compute the ground state of 
the Hamiltonian instead. 

This difficulty has been successfully overcome (see [9] 
and references therein) by redefining the reduced density 
matrix of the left block £ as: 

Pla' = EE^'*^(0o,^$(0c'./3, (3) 

pen I 

where a and a' label states in the left block C, /3 those of 
the right block TZ, and {$(/)}; is a set of, as of yet, unspec- 
ified states of the superblock CUTZ. The states $(/) are 
said to be targetted by the DMRG algorithm. Because of 
their inclusion in the reduced density matrix, these states 
will be obtained with a precision that scales similarly to 
the precision of the ground state in the static formulation 
of the DMRG. The relevance of the weights w/ appearing 
in Eq. (3) will be discussed in section II C. 



Which are the states <i>(Z) that need to be included in 
Eq. (3) to compute Eq. (1)? and \(j)2{t)) are cer- 

tainly needed. Since observables that include the ground 
state are ubiquitous, the ground state of the Hamiltonian 
needs to be targetted as well in most cases. But additional 
states need to be included in order to evolve to larger times, 
as we will now explain. 

Our implementation follows the time-step targetting 
procedure of Ref. [14]. We now introduce a small time 
T such that, for all t < t, Eq. (2) is accurate in the sense 
defined in, e.g., Ref. [22]. We consider a set of n„ times 
{tx}x, X = 0, 1, • • • rit, - 1, such that < t^+i, io = 0, and 
tn^-i = T. For simplicity, we assume from now on that 
= 102) = W) in Eq. (1). The state |0) is defined by the 
particular physics problem under investigation and we will 
consider particular examples in section III. States \4>{tx)) 
for each Q < x < can be obtained accurately from \4>) 
by using Eq. (2) since < r. To compute Eq. (1) for all 
t < T, we target the states \4>{tx)) and the ground state 
1^) as well. 

At this point it is instructive to consider a concrete class 
of states In a large class of problems these states are 
related to the ground state |^) of the Hamiltonian by 



3 - 4 - 5 - 6 
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for b local operators i?o,7r'(o)j • • • , act- 

ing on sites 7r'(0), 7r'(l), • • • ,7r'(&— 1) of the superblock, 
where 7r'(i) denotes the lattice site on which the operator 
-Bi,7r'(i) acts, as explained below. The extra index i indi- 
cates that the operators B can be different on different 
sites. Physical examples of the operators B will be given 
in section III. 

The sites 7r'(0), 7r'(l), • • • , 7r'(6 - 1) in Eq. (4) (as weh 
as sites 7r(0), 7r(l), • • • , 7r(a— 1) in Eq. (1)) will be consid- 
ered ordered in the way in which they appear as central 
sites [9] for the DMRG finite algorithm, as illustrated in 
Fig. 1. At a given stage of the computational procedure, 
if the central site of the DMRG algorithm is 7r'(0), then 
l^partiai^ = -Bo,-n-'(o) 1^) Can bc obtained. Next, we pro- 
ceed to the following site, and so on, until we reach site 
7r'(l), and apply i.e., I^ii'o'"'') = 

eventually reaching site ir'ib — 1), to complete the compu- 
tation of given by Eq. (4). Since in cases of physical 
interest the operators B are either bosons or fermions, a 
reordering is always possible due to their commutativity 
or anti-commutativity, yielding at most a minus sign. 

As the DMRG algorithm sweeps the entire lattice, the 
central sites change, leading to modified Hilbert spaces. 
Therefore, a procedure is required to "transport" the states 
\cf>) from one space to another. It is known [13] that the 
transformation needed to "transport" these states is the 
so-called wave-function transformation, that was proposed 



B4l¥> 

B,B^ |x|/> 
B3B,B^ |v|/> 



FIG. 1: Example of a state given by Eq. (4). The order in 
which the on-site operators are applied will depend on the order 
in which the sites appear as "central sites". In the example 
above, 6 = 3 and the operators are acting on sites 3, 4, and 
7. In the sweeping procedure shown, the order will be 7r'(0) = 
4,7r'(l) = 7,7r'(2) = 3, i.e., j0) = B3B7B4I1/'). 



by White [25] first in the context of providing a guess for 
the initial Lanczos vector to speed up the algorithm, but 
later found to be of applicability for other sub-algorithms. 

After the state |0) in Eq. (4) is computed, the DMRG 
algorithm operates for a few extra steps to better converge 
aU states li^it^)}, ^t^ < r. These states and aU DMRG 
transformations can be saved to disk, and later be used to 
compute the observables Eq. (1). 



C. Evolution to Arbitrary Times 

What happens to Eq. (1) for larger times, i.e., for times 
t > t7 Noting that |(/)(r -|- r)) = exp(?7?T)|0(r)) we can 
apply Eq. (2) to which in general reads [13]: 



E 

k,k',l,l',j 



(5) 

Then, we proceed by targetting the states {\4>itx + t))}^ 
for some time until they are converged. By applying this 
procedure recursively, we reach arbitrary times t as we 
sweep the finite lattice back and forth, and target {|(/'(ta; -I- 
t))}x in the general case. 

The speed of time advancement in the algorithm is con- 
trolled by two opposing computational constraints. If we 
advance times too fast by applying Eq. (5) too often, then 
convergence might not be achieved, or we might not have 
had the chance to visit all sites 7r(0), 7r(l), . . . , 7r(a — 1) 
to compute Eq. (1). Conversely, advancing too slowly 
would increase computational cost but produce no addi- 
tional data. Remember that when not advancing in time, 
states \(t){t + tx))^ are wave-function-transformed, as ex- 
plained in the previous section. 

We now explain the choice of the weights [14] that ap- 
pear in Eq. (3). Assume n to be the total number of states 
to be targetted, including the ground state. To give them 
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more prominence, we have chosen a weight of 217 for the 
ground state, and also for the $(/) vectors at the begin- 
ning and end of the r intervaL We have chosen a weight 
of fl for the rest. Then 217 x 3 + ri(n — 3) — 1 imphes 
= l/(n + 3). The algorithm does not appear much de- 
pendent on the choice of weights. However, irrespective of 
what the choice actually is, all mentioned states must have 
non-zero weights to avoid loss of precision for one or more 
states. 



D. Overview of the Implementation 

The DMRG+-^ code was introduced in Refs. [26, 27]. 
The extension of the code to handle the time evolution 
and computation of observables of the type represented by 
Eq. (1) was carried out with minimal refactoring. A Tar- 
getting interface was introduced, with two concrete classes, 
GroundStateTargetting, and TimeStepTargetting. The 
first handles the usual case, and is used even in the 
presence of time evolution during the so-called "infinite" 
DMRG algorithm, and during the finite algorithm before 
encountering the first site 7r'(0) in Eq. (4). 

A call to target . evolve (...) handles (i) the compu- 
tation of the vectors {|0(i + as needed, and (ii) their 
time evolution or, depending on the stage of the algorithm, 
their wave-function-transformation. When the target ob- 
ject belongs to the TimeStepTargetting class, the actual 
implementation of these tasks is performed by the mem- 
ber function evolve (. . .). When the target object is 
of class GroundStateTargetting the evolve (. . .) func- 
tion is empty. The call to this function is always done 
immediately after obtaining the ground state 1-0) for that 
particular step of the DMRG algorithm. 

File TimeStepTargetting. h is documented in place us- 
ing literate programming [28]. Further details about how 
to run the DMRG-I--I- code, and how to specify its input 
file are given in Appendix C. 



III. EXAMPLES 

A. One-site Excitations 

To test the accuracy of the time-dependent DMRG ap- 
proach explained in the previous section, we consider first 
the following problem. Consider the tight-binding model 
Hq = j ^tijcl^Cja, with tij a symmetric matrix, and 
with the observable we wish to calculate being 

X,,^{t) = (0|4e-*^*n,^e'^*c,^|v^), (6) 




time 



FIG. 2: Observable Xij-^{t), Eq. (6), for Hamiltonian Hq, with 
1 = 8 and j as indicated, on two geometries: (a) a 8x2 ladder, 
and (b) a 16-site chain. Circles and squares represent DMRG 
results, and solid lines exact results. Inset: Labelling of sites 
on a two-leg ladder. The highlighted site i = 8 is where the 
one-site excitation (creation of a "hole") was applied. Chain 
sites are labelled from left to right, starting at 0. 



where Ii/)) is the ground state of Hq. (We keep the usual 
notation tij for the matrix of hopping integrals in the con- 
text of a tight binding model [29], even though t is also 
used to denote time here.) 

This is equivalent to taking h — 1, B — c^, and 7r'(0) = i 
in Eq. (4); and a — 1, 7r(0) = j, and A — n^'wi Eq. (1). The 
physical interpretation for Xij^(t) is then clear: it provides 
the time-dependent expectation value of the charge density 
{nj^)(t) at site j over a state that, at time t = 0, is defined 
by creating a "hole-like" excitation in state \tjj) at site i. 
We assume that site i has been specified and is kept fixed 
throughtout this discussion. 

Xij {t) can be expressed in terms of the eigenvectors and 
eigenvalues of tij . For a half-filled lattice we have Xij (t) — 
R^Rj -\Tij (t) |2 , where R, and T,-j are given in Appendix A. 
Then, Xij{t) can be calculated numerically, and compared 
to DMRG results for this model on a 16-site chain, and on 
a 8x2 ladder, see Fig. 2. 

On the chain, the number of states kept for the DMRG 
algorithm was set to m = 200, which was found to give 
good accuracy for the ground state energy. On the ladder, 
m — 400 was used, which is a typical [30] m value to 
achieve good accuracy for the ground state energy, and 
static properties. For instance, in both cases, DMRG gives 
Xii{t = 0) = and Xij{t = 0) = 1/4 for i ^ j,as expected. 
As shown in Fig. 2, the use of these values of m yields an 
accurate time evolution. 
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B. Holon-Doublon 

Currently there is considerable interest in studying 
the feasibility of a new class of materials — the Mott 
insulators — for their possible use in photovoltaic devices 
and oxide-electronics in general. The crucial question un- 
der study is whether charge excitations in the Mott in- 
sulator will be able to properly transfer the charge into 
the metallic contacts, thus establishing a steady-state pho- 
tocurrent. Answering this question will require computa- 
tion of the out-of-equilibrium dynamics and the time evo- 
lution of the excitonic excitations produced by the absorp- 
tion of light by the material. 

The electron and hole created by light absorption are 
modeled by the state [31] 



,(1 - n,5)ct 



(7) 



where \ip) is the ground state, a and tr' are spin indices, 
and a = 1 — a denotes the spin opposite to a. A sum 
over (T and cr' is assumed in the equation above. This is 
equivalent to taking b = 2, Bq = c\n-,^ Bl = Ca-{1 — rig), 
7r'(0) — j, and 7r'(l) = i in Eq. (4). We assume that the 
sites i and j of the lattice have been specified and will 
remain fixed throughout this discussion. 

To model a Mott insulator we consider the Hubbard 
Hamiltonian [32-34] 



H = Hn 



(8) 



where the notation is as in Ref. [31], and we will drop the 
hat from the operators from now on. The hopping matrix 
t corresponds either to an open chain or to a two-leg ladder 
in the studies below. 



1. Density 

The time-dependent density at site p of state Eq. (7) is 



(9) 



which amounts to taking a = 1, A — n-^, and 7r(0) = p in 
Eq. (1). Consider Oj^i^p{t) = J2a Oj,i,p,a{t)- In the case of 
U = and half- filling we have (details are in Appendix B): 

ii p = i 

O,-,,p(0) = O,,,,p,t(0)-HO,,,,p,;(0) = <( 2 ifp = J (10) 

1 otherwise. 

The observable we test in this section is {'i e\np\'^ e) , 
which has a similar physical interpretation as Xij(t) 
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FIG. 3: Density at site p of the holon-doublon state Eq. (7) 
as a function of time, for (a) U — 0, and (b) U = 10. A 2x4 
ladder with = 1 and ty = 0.5 was used, containing 4 up and 



4 down electrons, 
i = 2, j = 4. 



The holon-doublon operator was applied at 



(Eq. (6)) in the holon-doublon case. Results for U = 
and U — 10 are shown in Fig. 3. At t = 0, the values 
of (*e|n-p|*e) given by Eq. (10) hold true in the J7 7^ 
case at half-filling. The time-evolution for interacting and 
non-interacting cases are, however, quite distinct, as in the 
case of the chain (see, e.g., Ref. [35]). 

Readers might want to know why we emphasize the non- 
interacting U = case. One obvious advantage of the 
U — case is that we can test the Krylov method, and 
indirectly the accuracy of the DMRG, against exact re- 
sults. In addition, the U term, at least when on-site, is 
not a major source of efficiency problems for the DMRG 
algorithm. 

To test our results for J7 ^ we have compared them 
to the Suzuki- Trotter method (not shown). We have also 
computed the small time expansion, and this is shown in 
Fig. 4. 



2. Double- occupation 
The double-occupation of state Eq. (7) is [31] 

iv.o...M). '*-"'?;ry"'i^> . (11) 



and amounts to taking a — 1, A — n^nj^, and 7r(0) = p in 
Eq. (1). 

Summarizing the operator equations obtained in sec- 
tion IIIBl, rii^Aij = 0, and fij^Aij = 0, where Aij 
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FIG. 4: Comparison of Oj.i,p{t) computed with DMRG (cir- 
cles and squares) and with a t — >■ expansion. For the 
latter Oj,i,p{t) = Oj,i,p(0) + aj^i^pt^ + 0{t'^), where aj.i.p = 
{cj>\HnpH — H^np\(j)). Same parameters as in Fig. 3b. 



is the operator defined in Eq. (Bl), n = 1 — n, and 
these equations also hold if we replace f by |. Then 
Nd{j,i,i,t = 0)(</)|(/)) = {AljUi^n^iAij) = 0, and 

DMRG results for [/ = and U — 10 are shown in 
Fig. 5. Also shown are exact results for U — 0. At i = 0, 
the double occupation at the doublon (p — j) and holon 
{p = i) sites are, respectively, Nd{j, i,p = j,t = 0) = 1 and 
^dij, i,p = i,t = 0) = for both interacting and noninter- 
acting cases. At the doublon site, the double occupation 
has a characteristic oscillating decay caused by the dy- 
namics of the holon-doublon pair within the system, also 
observed for the chain case [35]. 



IV. COMPUTATIONAL EFFICIENCY AND 
CONCURRENCY 

As in the static DMRG algorithm, the most computa- 
tionally intensive task of the time-step targetting DMRG 
algorithm is the computation of Hamiltonian connections 
between the system and environment blocks. The differ- 
ence is that now the lattice needs to be swept painstakenly 
to advance to larger and larger times. The scaling, how- 
ever, is linear with the number of finite sweeps, as long as 
the truncation m remains constant. 

This expensive task of building Hamiltonian connec- 
tions between system and environment blocks can be par- 
allelized [36-40]. Our implementation uses pthreads, a 
shared memory approach [45]. In percentage, the compu- 
tation speed-up is similar to the static DMRG case, and a 
discussion of the strong scaling can be found in Ref. [27] . In 
terms of wall-clock time, the speed-up is larger due to the 




FIG. 5: Double occupation at site p of the holon-doublon state 
Eq. (7) as a function of time, for (a) f/ = and (b) i7 = 10. A 
2x4 ladder with = 1 and ty = 0.5 was used, containing 4 up 
and 4 down electrons. The holon-doublon operator was applied 
at i = 2, j = 4. 



time-step targetting DMRG algorithm taking more time 
than the static version. 

The computation of target states could be parallelized 
easily, but whether serial or parallel, it is too fast to have 
substantial impact on the CPU times of production runs. 

The measurement of observables is a different matter. 
DMRGH — h computes observables post-processing, i.e., the 
main code saves all DMRG transformations, permutations, 
and quantum states to disk, and a second observer code 
reads the data from disk and computes observables as 
needed. We argue that post-processing is more advanta- 
geous than in-situ processing, whether for the static or for 
the time-dependent DMRG algorithm. 

First, a single run of the main code enables computa- 
tion of all observables. If, instead, one needed to make 
a decision on what observables to compute when running 
the main code, one would risk computing too much or too 
little. In the former case, computational resources and 
wall clock time would be wasted. In the latter case, the 
main run would have to be repeated, leading to vast redun- 
duncies because the observations are not computationally 
intensive compared to the main code. 

Moreover, by computing observables post-processing, we 
decouple the code, and enable scalable parallel compu- 
tations. For example, one-point observables of the form 
{(l)i{t)\Ai\(j)2{t)) are parallelized over i, and two-point cor- 
relations, such as {(t)i{t)\AiBj\4)2{t)), are parallelized over 
j, and cached over i. If is the number of sites of the 
lattice, the parallelization scales linearly up to almost N; 
the scaling is good but not perfect due to initialization 
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costs [41]. 

V. SUMMARY 

This paper explained in detail the implementation of 
the Krylov method for the real time evolution within the 
DMRG algorithm, using time-step targetting [13, 14]. We 
applied the method to a simple case of one-site excitations 
and found the method to be accurate. For the case of 
the holon-doublon excitation, we have extended to two- 
leg ladders the previous results obtained in chains. Our 
analysis has shown that the method is accurate as long as 
the underlying DMRG algorithm is accurate. Since Mott 
insulators are under study for its possible applicability to 
solar cells, the present results pave the way for their con- 
tinued study, now on more complex (but still quasi-one 
dimensional) geometries, such as ladders. 

We described computational tricks that can help de- 
crease the runtime. For example, we mentioned that 
shared memory parallelization with a few CPU cores can 
cut times by a factor of 2 or more. Parallelization works 
in the same way for time-dependent DMRG as it does for 
static DMRG, but helps more in the former case, due to 
runs taking longer. We also argued in favor of the post- 
processing of observables to speed-up production runs, and 
increase computational efficiency. 

Our implementation, DMRG-I--I-, is free and open 
source. It emphasizes generic programming using C-l — h 
templates, a friendly user-interface, and as few software 
dependencies as possible. DMRG-t-+ makes writing new 
models and geometries easy and fast, by using a generic 
DMRG engine. 

Appendix A: One-site Excitation in the 
Non-Interacting Case 

Let U be the matrix that diagonalizes Hq, so that c|.j. = 

\ i^'^A-T' ^^'^ '"-'^.t diagonal operators. Let E\ be 
the eigenvalues of Hq. 

For the rest of this appendix we omit After some 
algebra, and omitting the sums over duplicated indices: 

X^J{t) = U*^Ul^U,,^,U,,y{ule-''''ulu^e'"'ux'). (Al) 

The ground state of Hq is made up of iV^ filled levels, 
up to the Fermi energy, and particle-hole excitations are 
eigenstates of H (or, conversely, the excited states of Hq 
are particle- hole excitations). Then, the A'— th level of \(j>) 
is occupied and u\'\(j)) is an eigenstate of H with a hole at 
A'. Applying this reasoning multiple times, the final result 

is X,j{t) = R,Rj - |^y(^)|^ where R, = J^'x \U^,x\^, and 



Tij{t) — \Uj,xexp{~iE\t), where the prime over 

the sumation means sum only over occupied states A. 



Appendix B: Holon-Doublon for a Non-Interacting 
System 

The goal of this appendix is to compute Eq. (9) when 
U — and i j- Let N be the number of sites, let N be 
even, and let there be = iV^ = N/2 electrons. Then 
(m^) = (n^) = 5, Vi. 

Let Mij^ = (ni-fTij^) be an equation between real num- 
bers, and Hi = rii-^ + rij^ be an operator equation. From 
{rii^) — {nil) — 5' follows that 2Nij^ + 5 ~ (nirij). 
Another straightforward result is Oj,i,i,^(0) = 0. 

For an operator A let us define 'D{A) = (A^A), and let 

Aij = Cicr(l - n^a)clj„,nj-,, (Bl) 
then after some algebra: 

V{A,,) EE V,, = ^ - 2M,^ + 4A/-2^ + 2|A',,^|^ (B2) 

where A'^^ = (cl^Cj^)- 

By writing = 1 - Cp^c^^, one gets Oij j,t(0) = 1 ~ 
term. It is straightforward to prove that the second term 
vanishes and thus: Oj^ij^^{0) = 1. 

Now consider p ^ i and p ^ j. If we expand in the basis 
that diagonalizes the Hamiltonian we arrive to: 

O,„,p,t(0) = M +Y,{4JA^M{n\np^)l'D^„ (B3) 

n>0 

where the sum in n is over excited states. The second term 
is non-negative. This follows because A\-Aij and rip^ com- 
mute, admiting a common basis where both are diagonal. 
Moreover, all eigenvalues of AljAij are non-negative, as 
are all those of Up^. Then, the second term in Eq. (B3) is 
non-negative. 

Also, the first term of Eq. (B3) is i. Then Vp ^ i,j 

Oj,i,p,-f (0) = 5 + ^ijp^ where r^p > 0. 
If we sum over all sites, 

(4^,, ^ nr^)/V,, = {Al^A,N^)/n,, - ^. (B4) 

r 

There are N — 2 sites p such that p ^ Putting it all 
together we get: 

2^ Oj,i,r.t{0) = 1 + + 2^ = ^ + 2^ njp, 

(B5) 
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which, accordmg to Eq. (B4) has to be equal to It fol- 
lows that J2p^i,j "^ijp = 0' implying Vijp = Q\/p^i,j since 
we knew the rs were non- negative. Hence Oj^i^p^'^{Q) = ^, 
Vp ^ i, j 

Appendix C: The DMRG++ Code 

The required software to build DMRG++ is: (i) GNU 
C++, and (ii) the LAPACK and BLAS libraries [42]. 
These libraries are available for most platforms. The con- 
figure. pi script will ask for the LDFLAGS variable to pass to 
the compiler/linker. If the Linux platform was chosen the 
default /suggested LDFLAGS will include -llapack. If the 
OSX platform was chosen the default/suggested LDFLAGS 
will include -framework Accelerate. For other platforms 
the appropriate linker flags must be given. More informa- 
tion on LAPACK is here: http://netlib.org/lapack/. 

Optionally, make or gmake is needed to use the Makefile, 
and perl is only needed to run the configure.pl script. 

To build and run DMRG++: 

cd src 

perl conf igure . pi 

(please answer questions regarding model, and 
choose TimeStepTargetting) 
make dmrg; make observe 



. / dmrg input . inp 
./observe input . inp time 

The perl script configure.pl will create the 
files main.cpp, Makefile and observe. cpp. Ex- 
ample input files for one-site excitations are in 
TestSuite/inputs/inputS . inp, and for holon-doublon 
excitations in TestSuite/inputs/inputlO . inp. These 
files can be modified and used as input to run the 
DMRG++ program. Further details can be found in the 
file README in the code. 

The FreeFermions code found at http://www.ornl. 
gov/'gEl/FreeFermions/ can be used to compute prop- 
erties of non- interactive systems. 
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